Computation and visualization of Casimir forces in arbitrary geometries: 
non-monotonic lateral-wall forces and failure of proximity- force approximations 
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We present a method of computing Casimir forces for arbitrary geometries, with any desired ac- 
curacy, that can directly exploit the efficiency of standard numerical-electromagnet ism techniques. 
Using the simplest possible finite- difference implementation of this approach, we obtain both agree- 
ment with past results for cylinder-plate geometries, and also present results for new geometries. 
In particular, we examine a piston-like problem involving two dielectric and metallic squares sliding 
between two metallic walls, in two and three dimensions, respectively, and demonstrate non-additive 
and non-monotonic changes in the force due to these lateral walls. 



PACS numbers: 

Casimir forces arise between macroscopic objects due 
to changes in the zero-point energy associated with quan- 
tum fluctuations of the electromagnetic field [1 j. This 
spectacular effect has been subject to many experimental 
validations, as reviewed in Ref. 2, All of the experiments 
reported so far have been based on simple geometries 
(parallel plates, crossed cylinders, or spheres and plates). 
For more complex geometries, calculations become ex- 
tremely cumbersome and often require drastic approx- 
imations, a limitation that has hampered experimental 
and theoretical work beyond the standard geometries. 

In this letter, we present a method to compute Casimir 
forces in arbitrary geometries and materials, with no un- 
controlled approximations, that can exploit the efficient 
solution of well-studied problems in classical computa- 
tional electromagnetism. Using this method, which we 
first test for geometries with known solutions, we predict 
a non-monotonic change in the force arising from lateral 
side walls in a less-familiar piston-like geometry (Fig.j2|. 
Such a lateral-wall force cannot be predicted by "addi- 
tive" methods based on proximity- force or other purely 
two-body-interaction approximations, due to symmetry, 
and it is difficult to find a simple correction to give a 
non-monotonic force. We are able to compute forces for 
both perfect metals and arbitrary dispersive dielectrics, 
and we also obtain a visual map of the stress tensor that 
directly depicts the interaction forces between objects. 

The Casimir force was originally predicted for parallel 
metal plates, and the theory was subsequently extended 
to straighforward formulas for any planar-multilayer di- 
electric distribution s(x,uj) via the generalized Lifshitz 
formula [3]. In order to handle more arbitrary geome- 
tries, two avenues have been pursued. First, one can em- 
ploy approximations derived from limits such as that of 
parallel plates; these methods include the proximity- force 
approximation (PFA) and its refinements [4], renormal- 
ized Casimir- Polder [5] or semi-classical interactions [6], 
multiple-scattering expansions [7 , classical ray optics [8] , 



and various perturbative techniques [9| HOj . Such meth- 
ods, however, involve uncontrolled approximations when 
applied to arbitrary geometries outside their range of ap- 
plicability, and have even been observed to give quali- 
tatively incorrect results [TTJ[T2]. Therefore, researchers 
have instead sought numerical methods applicable to ar- 
bitrary geometries that converge to the exact result given 
sufficient computational resources. One such method 
uses a path-integral representation for the effective ac- 
tion [13 , and has predicted the force between a cylin- 
der and a plane or between corrugated surfaces. Ref. [13] 
uses a surface parameterization of the fields coupled via 
vacuum Green's functions, requiring 0(N 2 ) storage and 
0(N 3 ) time for N degrees of freedom, making scaling 
to three dimensions (3d) problematic. Another exact 
method is the "world-line approach" [I2j [14] [15] , based 
on Monte-Carlo path-integral calculations. (The scaling 
of the world-line method involves a statistical analysis, 
determined by the relative feature sizes in the geome- 
try, that is beyond the scope of this Letter.) Further- 
more, the methods of Ref. IT3l and Ref. [Ml have currently 
only been demonstrated for perfect-metallic z-invariant 
structures — in this case, the vector unknowns can be de- 
composed into TE (E • z = 0) and TM (H • z = 0) scalar 
fields with Dirichlet or Neumann boundary conditions — 
although generalizations have been proposed [16] . Here, 
we propose a method based on evaluation of the mean 
stress tensor via the fluctuation-dissipation theorem, 
which only involves repeated evaluation of the electro- 
magnetic imaginary- frequency Green's function. For a 
volume discretization with N degrees of freedom and 
an efficient iterative solver, this requires O(N) stor- 
age and 0(N 2 ~ 1 ^ d ) time in d dimensions. Further- 
more, because evaluation of the Green's function is such 
a standard problem in classical computational electro- 
magnetism, it will be possible to exploit many develop- 
ments in fast solvers, such as finite-element, or boundary- 
element methods [T7] , To illustrate the method, our 
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initial implementation is based on the simplest-possible 
finite-difference frequency-domain (FDFD) method, as 
described below. 

As derived by Dzyaloshinskh et al. [I], the net Casimir 
force on a body can be expressed as an integral over any 
closed surface around the body of the mean electromag- 
netic stress tensor (T^), integrated over imaginary fre- 
quencies ou = iw: 

Fi= r**E K Y,(T. 3 (r,iw))dS 3 . (1) 

For a 3d z-invariant structure, the z integral is replaced 
by an integral over the corresponding wavevector, result- 
ing in a net force per unit length. The stress tensor is 
defined as usual by: 

(T ZJ (r,iw)) = (Hi(r) Hj(r)) - % ]T (H k (r) H k (r)) 



+ e(r, iw) 



(E^E^-h^iE^EUr)) 



(2) 



The connection to quantum mechanics arises from the 
correlation functions of the fluctuating fields, such as 
(EiEj), given via the fluctuation-dissipation theorem in 
terms of the imaginary-cj Green's function Gij(iw; r — r'): 

(E i (r)E j (r')) = V ^G ij (iw;r-r') (3) 
(fl i (r)JJ i (r / )) = -(Vx) l( (V'x) jm G, m (w;r-r'), (4) 
where the Green's function G^ solves the equation: 



V x V x 



-e(r, iw) 



Gj(iwir-r') = e^r-r') (5) 



for a unit vector ej in the j direction, and obeys the usual 
boundary conditions on the electric field from classical 
electromagnetism. (The above expressions are at zero 
temperature; the nonzero-temperature force is found by 
changing J dw in Eq. [I] into a discrete summation pQ.) 
Although the Green's function (and thus T^ ) is formally 
infinite at r = r', this divergence is conventionally re- 
moved by subtracting the vacuum Green's function; in 
a numerical method with discretized space, as below, 
there is no divergence and no additional regularization 
is required. (The vacuum Green's function gives zero net 
contribution to the dS integral, and therefore need not 
be removed as long as the integrand is finite.) 

Historically, this stress-tensor expression was used to 
derive the standard Lifshitz formula for parallel plates, 
where Gij is known analytically. However, it also forms 
an ideal starting point for a computational method, be- 
cause the Green's function for arbitrary geometries is 
routinely computed numerically by a variety of tech- 
niques [17]. Furthermore, the problem actually becomes 



easier for an imaginary uj. First, for an imaginary uj, the 
linear operator in Eq. [5] is real-symmetric and positive- 
definite for w 7^ 0, since the dielectric function e(uj) is 
purely real and positive along the imaginary-o; axis for 
physical materials without gain, due to causality. Second, 
the imaginary-o; Green's function is exponentially decay- 
ing rather than oscillating, leading to a well-behaved non- 
oscillatory integrand in Eq. [T] 

To illustrate this method, we employed the simplest 
possible computational technique: we perform a FDFD 
discretization of Eq. [5] with a staggered Yee grid [T8] 
and periodic boundaries, inverting the linear operator by 
a conjugate-gradient method. The presence of discon- 
tinuous material interfaces degrades second-order finite- 
difference methods to only first-order accuracy, and the 
uniform spatial resolution is also suboptimal, but we 
found FDFD to be nevertheless adequate for small 2d 
geometries. The periodicity leads to artificial "wrap- 
around" forces that decay rapidly with cell size L (at least 
as 1/L 3 in 2d and 1/L 4 in 3d); we chose cell sizes large 
enough to make these contributions negligible (< 1%). 

The computational process is as follows: pick some sur- 
face/contour around a given body, evaluate the Green's 
function for every grid point on this surface in order to 
compute the surface integral of the stress tensor, which 
is then integrated over w by adaptive quadrature. 




FIG. 1: Casimir force between a radius- R cylinder and a 
plate (inset), relative to the proximity-force approximation 
i^PFA, vs. normalized separation a/R. The solid lines are 
the Casimir force computed in Ref. 19 for TE (gray) and 
TM (blue) polarizations, along with results computed by our 
method with a simple finite- difference discretization (gray 
squares). Error bars were estimated for some data points 
by using computations at multiple spatial resolutions. Inset 
shows interaction stress tensor A(T XX ) at a typical imaginary 
frequency w = 2-kc/cl, where red indicates attractive stress. 

Before we attempt to study new geometries with our 
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method, it is important to check it against known re- 
sults. The simplest cases, of parallel metallic or dielec- 
tric plates, of course match the known result from the 
Liftshitz formula and are not reproduced here. A more 
complicated geometry, consisting of a perfect metallic 
cylinder adjacent to a perfect metallic plate in 3d, was 
solved numerically by Ref. [19l to which our results are 
compared in Fig. [T] Ref. [19] used a specialized Fourier- 
Bessel basis specific to this cylindrical geometry, which 
should have exponential (spectral) convergence. Our use 
of a simple uniform grid was necessarily much less effi- 
cient, especially with the first-order accuracy, but was 
able to match the Ref. [191 results within ~ 3% using rea- 
sonable computational resources. A simple grid has the 
advantage of being very general, as illustrated below, but 
other general bases with much greater efficiency are pos- 
sible using finite-element or boundary-element methods; 
the latter, in particular, could use a spectral Fourier ba- 
sis analogous to Ref. [19] and exploit a fast-multipole or 
similar 0(N log N) solver technique [T7] . 

Also shown, in the inset of Fig. [I] is a plot of the 
interaction stress-tensor component A(T XX ) at a typical 
imaginary frequency w = 2irc/a. By "interaction" stress- 
tensor A(T^), we mean the total (T^ ) of the full geome- 
try minus the sum of the (T^-)'s computed for each body 
in isolation. Here, the stress tensors of the isolated cylin- 
der and plate have been subtracted, giving us a way to 
visualize the force due to the interaction. As described 
below, such stress plots reveal the regions in which two 
objects most strongly affect one another, and therefore 
reveal where a change of the geometry would have the 
most impact. (In contrast, Ref. [12] plots an interaction- 
energy density that does not directly reveal the force, 
since the force requires the energy to be differentiated 
with respect to a. For example, Ref. [12js subtracted en- 
ergy density apparently goes nearly to zero as a metallic 
surface is approached, whereas the stress tensor cannot 
since the stress integration surface is arbitrary.) 

We now consider a more complicated geometry in 
which there are interactions between multiple bodies: a 
3d "piston" -like structure, shown in the inset of Fig. |2j 
consisting of two z-invariant metal s x s squares sepa- 
rated by a distance a from one another (here, s = a) 
and separated by a distance ft from infinite metal plates 
on either side. We then compute the Casimir force per 
unit z between the two squares as a function of the sep- 
aration ft. The result for perfect conductors is shown 
in Fig. |2j plotted for the TE and TM polarizations and 
also showing the total force. (Error bars are not shown 
because the estimated error is < 1%.) In the limit of 
ft — > 0, this structure approaches the "Casimir piston," 
which has been solved analytically [20j [21] (and also in 
2d for the TM polarization [22]). Our results, extrapo- 
lated to ft = 0, agree with these results to within w 2% 
(although we have computational difficulties for small ft 
due to the high resolution required to resolve a small 
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FIG. 2: Casimir force per unit length between metal squares 
F I Fpfa 5 vs. distance from metal plate ft (inset), normalized 
by the total TE+TM force per unit length obtained using the 
PFA, Fpfa = hcs((3)/4807ra 4 . The total force is plotted 
(black squares) along with the TE (red dots) and TM (blue 
circles) contributions. 



feature in FDFD). For ft > 0, however, the result is sur- 
prising in at least two ways. First, the total force is 
non-monotonic in ft, due to a competition between the 
TE and TM contributions to the forces. Second, the ft 
dependence of the force is a lateral effect of the parallel 
plates on the squares, which would be zero by symmetry 
in PFA or any other two-body-interaction approxima- 
tion. Although lateral-wall effects can clearly arise qual- 
itatively in various approximations, such as in ray optics 
or in PFA restricted to "line-of-sight" interactions, non- 
monotonicity is more surprising [24 . Also, in the large-ft 
limit, the force remains different from PFA due to finite- s 
"edge" effects [12], which are captured by our method. 

Our method is also capable, without modification, of 
handling dielectric materials. This is demonstrated in 
Fig. [3j where the Casimir force is shown for the case 
where the squares are made of gold, using the experimen- 
tal Drude e(uo) from Ref. 23 for a separation a = l/am. 
Here, our calculation is for a purely 2d geometry (equiv- 
alently, 3d restricted to z-invariant fields/currents), and 
for comparison we also plot the corresponding 2d force 
for perfect-metal squares (although the two cases are 
normalized differently as described in the caption). As 
might be expected, the dielectric squares have a weaker 
interaction than the perfect-metal squares, but are still 
non-monotonic. Note also the qualitative similarity be- 
tween the perfect-metal results of Figs. [2] and [3} reflecting 
the fact that the force contributions are dominated by 
the zero-wavevector (z-invariant) fields. Extrapolated to 
ft = 0, the perfect-metal TM force agrees with the known 
analytical result [22] to within « 3%. 




FIG. 3: Solid lines: Casimir force between 2d gold squares 
F j Fpfa i vs. distance from metal plate h (inset), using ex- 
perimental s{uj) [23] , normalized by the total force obtained 
using the PFA. (Here, the PFA force is computed for x-infinite 
gold slabs). The total force is plotted (black squares) along 
with the TE (red dots) and TM (blue circles) contributions. 
Dashed lines: force for 2d perfect-metal squares, normalized 
by the perfect-metal PFA force Fpfa = hcs ( (3) /8tv a 3 . 
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FIG. 4: (a-f): TM stress map of the 2d-analogous geome- 
try of Fig. [2] for various h. The intearaction stress tensors 
(T xx ) (left) and (T xy ) (right) for: (a),(d): h = 0.5a; (b),(e): 
h = a; and (c),(f): h = 2a, where blue/white/red = repul- 
sive / zero / attractive. 



To further explore the source of the ^-dependence, we 
plot the TM interaction-stress maps A(T XX ) and A(T xy ) 
in Fig. [4j for the 2d perfect-metal squares from Fig. [3] 
The stress plots of Fig. [4] are computed at a typical fre- 
quency w = 27rc/a, and for varying distances from the 
metal plates (h = 0.5, 1.0, 2.0). As shown, the mag- 
nitudes of both the xx (a-c) and xy (d-f) components 
of the stress tensor change dramatically as the metal 
plates are brought closer to the squares. For example, 
one change in the force integral comes from T xy , which 
for isolated squares has an asymmetric pattern at the 
four corners that will contribute to the attractive force, 
whereas the presence of the plates induces a more sym- 
metric pattern of stresses at the four corners that will 
have nearly zero integral. This results in a decreasing 
TM force with decreasing h. Because stress maps indi- 
cate where bodies interact and with what signs, it may 
be useful in future work to explore whether they can be 
used to design unusual behaviors such as non- additive, 
non-monotonic, or even repulsive forces. 
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